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ABSTRACT 

We  have  examined  the  effect  of  seismic  anisotropy  on  the  scattering  of  surface  waves  in 
actively-deforming  continental  regions,  many  of  which  occur  in  areas  of  nonproliferation  con¬ 
cern.  We  have  approached  the  problem  from  three  angles,  1)  we  have  developed  waveform 
inversion  methods  for  Love- Rayleigh  coupling  in  seismic  data  from  the  Tibetan  Plateau  and 
Tien  Shan  region,  2)  we  have  developed  a  flat-layered  surface  wave  code  to  study  the  effect  of 
crust  and  upper-mantle  anisotropy  on  crustal  resonances,  and  3)  we  have  advanced  wavelet- 
based  signal  processing  methods  to  estimate  the  polarization  of  coherent  seismic  energy 
across  a  broadband  seismic  array.  In  data  from  the  recent  portable  broadband  PASSCAL 
deployment  in  Tibet,  we  found  clear  evidence  of  long-period  Love-to- Rayleigh  scattering  that 
is  best  explained  by  lateral  gradients  of  anisotropy  in  the  upper  mantle,  developed  in  the 
course  of  the  compression  and  uplift  of  the  plateau.  We  developed  a  modal  summation  tech¬ 
nique  for  waveform  perturbations  that  allowed  us  to  infer  that  both  S  and  P  wave  anisotropy, 
with  a  ratio  consistent  with  the  mineral  alignment  of  sheared  peridotite,  in  the  depth  range 
100-300  km  are  necessary  to  fit  the  observed  data.  Using  a  plane-layered  geometry,  we  are 
investigating  the  effect  of  crust  and  upper  mantle  anisotropy  on  Love-Rayleigh  coupling  in 
the  5-30  second  period  range.  Most  previous  studies  have  posited  a  horizontal  or  vertical  axis 
of  symmetry  for  seismic  anisotropy,  but  our  code  can  incorporate  a  tilted  axis  of  symmetry, 
which  enhances  coupling  and  scattering  effects.  Such  dipping  geometries  could  develop  in 
the  crust  as  a  result  of  thrust  ramps  in  compressive  tectonic  regimes,  either  in  the  form  of 
oriented  crystals  in  shear  zones  or  fine-layering  of  isotropic  material.  We  have  developed 
a  multiple-station  wavelet-based  method  for  detecting  frequency- dependent  polarization  in 
noisy  seismic  data.  The  wavelet  method,  unlike  the  moving- window  Fourier  transform,  scales 
the  time  interval  analysed  to  the  period  of  the  presumed  signal.  We  have  applied  this  method 
successfully  thusfar  to  Terrascope  data  from  Southern  California,  and  have  documented  a 
strong  refraction  of  intermediate-period  surface  waves  along  the  ”big-bend”  portion  of  the 
San  Andreas  Fault. 
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OBJECTIVE 


We  wish  to  determine  the  effect  of  seismic  anisotropy,  within  the  crust  and  uppermost 
mantle,  on  surface  waves  with  periods  between  10  and  30  seconds.  In  particular,  we  wish 
to  characterize  the  scattering  caused  by  lateral  variations  in  anisotropy  in  regions  of  past 
and  present  lithospheric  deformation.  This  path- dependent  scattering  can  distort  seismic 
waveforms,  and  make  more  difficult  the  discrimination  of  explosions  from  earthquakes.  We 
plan  to  identify  regions  of  the  Earth  where  scattering  associated  with  lithospheric  anisotropy 
is  important.  Data  from  regional  networks  and  arrays  located  in  Eurasia  and  the  US  will 
also  be  used  in  experiments  to  extract  ’scattered’  waveforms  from  the  principal  Love  and 
Rayleigh  wavepackets.  In  addition  to  standard  stacking/beamforming  techniques,  we  will 
develop  a  frequency-dependent  variant  of  the  MUSIC  algorithm  to  identify  coherent  surface 
waves  that  traverse  an  array.  We  plan  to  test,  in  selected  areas,  whether  anisotropy  inferred 
from  the  ’splitting’  of  1-20  second  shear  waves  is  sufficient  to  explain  observed  Love  to 
Rayleigh  scattering. 

TECHNICAL  PROPOSAL:  RESEARCH  ACCOMPLISHED 

We  have  developed  codes  for  the  dispersion  characteristics  of  crustal  surface  waves  in  a 
1-D  anisotropic  media,  with  anisotropy  that  is  hexagonally  symmetric,  but  with  an  axis  of 
symmetry  that  is  arbitrary.  Previous  investigators  have  restricted  attention  to  symmetry 
axes  that  are  either  vertical  or  horizontal,  but  a  tilted  axis  of  symmetry  is  plausible  in 
the  uppermost  crust  and  perhaps  near  the  Moho,  where  thrusting  and  extension  involves 
a  combination  of  vertical  and  horizontal  deformation.  A  modification  of  this  code  can  be 
used  to  examine  the  behavior  of  teleseismic  body  waves  that  propagate  upward  through  an 
anisotropic  crustal  structure.  Enhanced  coupling  between  P-SV  and  the  SH  particle  motion 
occurs  for  waves  that  have  significantly  non-vertical  incidence,  in  particular  between  P  and 
SH  motion. 

We  would  like  to  express  the  interaction  of  upgoing  and  downgoing  plane  waves  in  a 
weakly  anisotropic  layered  medium  with  hexagonal  symmetry.  Hexagonal  symmetry  implies 
a  stress-strain  relation  that  possesses  an  axis  of  symmetry  w.  Elastic  properties  are  constant 
within  each  of  a  stack  of  layers  over  a  uniform  halfspace,  but  the  axis  of  symmetry  can.  vary 
among  the  layers.  We  express  the  elastic  properties  as  a  function  of  depth  as  A(r),  where  Aijki 
is  the  fourth-order  stress-strain  tensor.  If  the  axis  of  symmetry  is  horizontal,  we  can  express 
the  azimuthal  dependence  of  the  squared  P  and  SV  velocities  for  horizontal  propagation 
in  terms  of  the  angle  £  from  w,  according  to  formulas  similar  to  the  head- wave  formulas  of 
Backus ,  [1965],  Crampin  [1977]  and  Park  [1993]: 

pa2(()  =  A  -f  B  cos  2£  +  C  cos  4£ 
pj3\i)  =  D  E  cos  2£. 

If  density  perturbations  are  neglected,  knowledge  of  A,  B,  C,  D,  E  is  sufficient  to  determine 
the  stress-strain  tensor  [Shearer  and  Orcutt ,  1986].  Park  [1993]  showed  how  these  azimuthal 
relations  generalize  to  other  orientations  of  w.  We  assume  a  fiat  earth,  z  =  0  at  the  free 
surface,  and  z  increasing  downward.  We  assume  a  plane- wave  solution  of  the  form  U(x,  t)  = 
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u(x)e*(k'x  "*).  The  strain  tensor 

e  =  ^(Vu  +  (Vu)T)  =  ^(k  ®  u  -f  u  ®  k)e,^k'x_t*'^ 

In  each  layer,  the  elastic  tensor  can  be  expressed 

A  =  A&a  4  BAb  4-  C Ac  +  DA.d  +  EAe,  (2) 

where 

A^  =  1 0 1 

Ab  =  W0I  +  I0W 

Ac  =  8W  ®  W  —  I  ®  I  (3) 

Ad  =  (i5)I  ®  I  +  (/-<)1 0 1  —  21  ®  I 

A  e  =  2[(i5)AB  +  ( u)Ab  —  2Ab]  4-  Ad, 

where  £®’  is  the  tensor  product  operation,  and  W  =  w  ®  w  —  |I,  where  I  is  the  identity 
tensor.  The  permutation  (ij)  indicates  the  interchange  of  the  ith  and  jth  tensor  index  e.g., 
{(13)1  ®  =  SkjSn.  An  isotropic  elastic  tensor  A^°l  contains  only  terms  proportional  to 

the  isotropic  tensors  Ax  and  Ac,  as  neither  depends  on  w. 

Because  P,  SV  and  SH  motion  can  couple  at  all  interfaces  in  a  layered  anisotropic  structure 
with  general  axis  of  symmetry  w,  all  interacting  upgoing  and  downgoing  plane  waves  share 
the  same  horizontal  wavenumber.  Let  k  =  px  -f-  uz,  where  p  and  v  are  the  horizontal  and 
vertical  slownesses,  respectively.  In  a  general  anisotropic  medium,  there  are  three  upgoing 
and  three  downgoing  plane  waves.  For  fixed  p,  the  Christoffel  relations  lead  to  six  possible 
solutions  for  v,  as  determined  by  the  quadratic  eigenvalue  problem 

(p2T  -  pi  +  vpS  +  v2K)  •  u  =  0  (4) 

where 

T  =  (A-  B  +  C  -  D  +  E)x®x  +  {D  +  (2w2  -  1)E)I  +  (8  w2C  +  2E)w  ®  w 
-j-(B  —  4 C  —  2E)wx(v/  ®  x  +  x  ®  w) 

S  =  [A  —  B  +  C  —  D  +  E)(-k  ®  z  +  z  ®  x)  +  AwxwzEI  +  16wxwzCw  ®  w 

+(B  -4 C  -  2E)[wz(v?  ®  x  +  x  ®  w)  +  wx( w  ®  z  +  z  0  w)]  (5) 

R  =  (A- B  +  C-D  +  E)i®z  +  (D  +  {2w\  -  l)E)I  +  (8  w\C  +  2E)w  0  w 
+(B  —  AC  —  2E)wz(vf  ®  z  +  z  0  w) 


To  transform  (8)  into  a  linear  6x6  eigenproblem,  define  T  =  R  1-(p2T—pT)  and  S  =  pR  1-S 
to  obtain  (T  4-  yS  +  u2T)  •  u  =  0.  This  3x3  system  can  be  solved  with  the  6x6  system 


(T  -  vl)  -  U  = 


-S  -T 

I  0 


—  v 


I  0 
0  I 


(6) 


Solving  (6)  will  obtain  three  upgoing  and  three  downgoing  waves,  aside  from  exceptional 
cases  e.g.,  if  p  equals  the  P  or  S  slowness  in  a  layer  and  the  associated  wave  travels  horizon¬ 
tally.  (In  this  case  the  matrix  system  is  formally  defective.)  To  solve  the  equation  of  motion 
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in  a  stack  of  constant-property  layers,  one  matches  boundary  conditions  for  traction  and 
displacement  continuity  for  these  six  waves  at  each  interface.  We  obtain  surface- wave  modes 
by  forcing  an  evanescent  condition  in  the  halfspace,  propagating  the  stress-displacement  vec¬ 
tor  to  the  free  surface,  and  enforcing  a  zero-traction  condition.  There  axe  several  schemes 
for  doing  this.  We  have  adapted  the  computational  algorithm  of  Chen  [1993],  which  forces 
all  evanescent  waves  to  ’’decay”  away  from  interfaces  to  decrease  roundoff  error.  At  a  fixed 
frequency,  the  algorithm  searches  the  allowable  interval  of  horizontal  slowness  for  roots  to 
the  secular  function  that  describes  free-surface  traction.  The  computational  formulas  are 
quite  similar  to  the  Rayleigh- wave  formulas  of  Chen  [1993],  except  that  the  reflection  and 
transmission  matrices  are  3x3,  rather  than  2x2.  The  numerical  application  is  compli¬ 
cated  by  the  need  to  separate  travelling  modes  that  are  closely  spaced  in  phase  velocity.  To 
keep  track  of  solutions,  we  use  the  group  velocity  of  a  mode  at  a  frequency  /  to  predict 
a  phase- velocity  ’’window”  in  which  a  solution  should  be  found  at  frequency  /  +  Sf.  This 
also  reduces  the  number  of  propagator-evaluations  of  the  surface  traction  condition,  each  of 
which  requires  solution  of  a  6  x  6  eigenvector  system  in  each  layer.  This  type  of  interpolator 
is  not  always  successful,  unfortunately,  because  the  coupling  of  Rayleigh  and  Love  motion 
can  cause  ’’bow-ties”  in  the  dispersion  curves  where  surface- wave  modes  mix  in  near-equal 
amounts  -  see  Figure  1. 

We  have  developed  and  coded  a  path-integral  formalism  to  invert  Love-to- Rayleigh  and 
Rayleigh-to-Love  scattered  waveforms.  Although  we  plan  to  apply  these  techniques  to  flat- 
layer  geometry  for  crustal/lithospheric  surface  waves,  our  earlier  work  made  natural  the 
initial  development  of  a  spherical-earth  formalism  using  coupled  free  oscillations.  A  problem 
when  calculating  coupled-mode  synthetics  in  a  3-D  earth  model  is  the  sum  over  all  singlets 
associated  with  spherical  haxmonic  multiplets.  By  using  the  summation  theorem  of  spher¬ 
ical  harmonics,  we  sum  the  singlets  analytically  and  calculate  the  synthetics  by  summing 
multiplets  rather  than  singlets.  We  expand  the  elastic  tensor  to  a  general  6x6  matrix  to 
handle  arbitrary  orientations  of  the  anisotropic  fast  axis.  This  technique  enables  us  not  only 
to  calculate  synthetics  quickly  in  a  3-D  anisotropic  structure  but  also  to  establish  an  in¬ 
verse  algorithm  because  of  the  linear  property  of  differential  seismograms.  These  synthetics 
compare  well  with  3-D  strong  Born  synthetics,  which  implies  that,  for  the  models  we  have 
considered,  Love-Rayleigh  scattering  occurs  primarily  along  the  great-circle  path.  We  use 
the  great-circle  technique  to  fit  waveforms  recorded  in  a  recent  broadband  PASSCAL  exper¬ 
iment  in  Tibet  (Figure  2).  We  find  that  the  observed  waveform  anomalies  in  the  southern 
part  of  the  Plateau  can  be  well  explained  with  a  strong  anisotropic  gradient  zone  beneath  the 
central  Plateau.  An  additional  anisotropic  zone  in  central  Alaska  can  be  used  to  fit  details  in 
the  Tibet  data  better.  The  presence  of  such  an  anisotropic  zone  is  supported  by  a  quasi-Love 
wave,  observed  at  College,  Alaska,  in  data  from  the  1/16/95  Kobe,  Japan  earthquake. 

We  have  also  extended  the  multiple  wavelet  polarization  method  to  apply  to  three- 
component  array  data.  We  have  been  testing  this  method  on  data  from  the  Terrascope 
broadband  regional  network  in  Southern  California.  The  wavelet  method  enables  the  user 
to  pick  apart  the  time-frequency  behavior  of  seismic  data,  and  the  SVD-polarization  code 
enables  the  user  to  correlate  motion  between  stations,  preserving  the  small  time  shifts  and 
amplitude/polarization  variability  that  would  be  averaged  out  in  simple  beamforming.  For 
surface  waves  in  the  40-100  sec  period  range,  there  is  considerable  variability  in  amplitude 
and  polarization  within  the  roughly  250  km  aperture  of  the  Terrascope  network.  There  is 
also  evidence  for  significant  Love-to-Rayleigh  forward  scattering  in  the  40-100  sec  period 
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Figure  1.  Dispersion  curves  for  hybrid  Love-Rayleigh  surface  waves  in  an  anisotropic 
crustal  waveguide.  The  velocity  model  is  a  modified  form  of  the  one- dimensional  velocity 
model  used  to  locate  earthquakes  in  Southern  California.  Anisotropy  with  a  horizontal  axis  of 
symmetry  at  45°  angle  to  the  propagation  is  assumed,  with  4%  peak-to-peak  variation  in  P- 
velocity  and  2%  peak-to-peak  variation  in  S-velocity,  both  with  sin  26  azimuthal  dependence. 
The  anisotropy  is  prescribed  only  in  a  5-km  thick  layer  at  the  base  of  the  crust,  consistent 
with  the  crustal  extension  common  in  the  western  US.  The  Love  and  Rayleigh  dispersion 
curves  tend  to  couple  strongly  where  they  cross  in  the  phase-velocity  graph,  making  ’’bridges” 
between  pairs  of  group  velocity  curves.  The  eigenfunctions  for  these  dispersion  branches  (not 
shown)  exhibit  significant  mixing  of  Love  and  Rayleigh  particle  motion,  even  far  from  these 
crossing  points. 
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LSA  Data:  6/28/92  Landers  event  Lithosphere  Anisotropy  (P+S) 


Figure  2.  The  upper  right  panel  shows  broadband  data  from  28  June  1992  Landers, 
California  event,  recorded  at  CDSN  station  LSA  (Lhasa,  Tibet).  The  lower  right  panel 
shows  the  coupled-mode-sum  synthetics  with  only  S  wave  (2%)  anisotropy,  with  an  east-west 
horizontal  symmetry  axis,  at  100-300  km  depth,  consistent  with  lithospheric  compression 
beneath  the  Tanggula  Shan  mountain  range  in  southern  Tibet.  The  lower  left  panel  shows 
the  synthetics  with  only  P  wave  (6%)  anisotropy.  The  upper  left  panel  shows  the  synthetics 
combined  with  6%  P  anisotropy  and  2%  S  anisotropy. 
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range,  which  indicates  lithospheric,  perhaps  crustal,  scattering.  Although  rough  isotropic 
lateral  structure  can  scatter  Love  to  Rayleigh  motion,  significant  scattering  usually  involves  a 
strong  deflection  of  the  wave.  Anisotropic  structure  is  a  better  model  for  forward  scattering, 
as  the  spatial  distortion  of  the  stress-strain  relation  better  couples  the  Love  and  Rayleigh 
particle- motion  polarizations. 

We  have  looked  at  seismic  data  recorded  by  Terrascope  from  Pacific  earthquakes  in  a 
broad  azimuthal  range,  to  examine  forward-scattering  across  the  array.  For  the  24  May 
1994  strike-slip  earthquake  in  Taiwan  (Ms= 6.6,  depth=33  km),  surface  waves  intersect  the 
California  coastline  at  a  low  angle  to  the  North- America/Pacific  plate  boundary,  and  so  are 
susceptible  to  severe  multipathing  and  scattering.  At  the  inland  station  GSC  (Goldstone, 
CA)  the  wavelet  transform  of  the  three  components  (vertical,  radial,  transverse)  reveals 
significant  radial-component  energy  for  the  Love  wave,  and  transverse-component  energy  for 
the  Rayleigh  wave.  A  clear  quasi-Love  scattered  wave  follows  the  Love  wave  closely  on  the 
vertical  component,  in  the  period  range  40-100  seconds.  Wavelet  polarization  analysis  of 
the  GSC  record  (Figure  3)  shows  that,  though  long-period  (>  100  sec)  Love  polarization 
is  undisturbed,  the  intermediate-period  polarization  suffers  refraction  and  some  elliptical 
motion  on  the  horizontals.  The  projection  of  correlated  motion  onto  the  radial- vertical  plane 
shows  low-amplitude  elliptical  motion,  consistent  with  the  quasi-Love  wave’s  ’’Rayleigh” 
nature.  The  first  pulse  of  Rayleigh  energy  is  strongly  refracted  at  intermediate  frequencies, 
by  roughly  45  degrees! 

We  have  solved  for  polarized  signals  that  are  correlated  network-wide  using  a  singular 
value  decomposition  (SVD)  of  a  matrix  of  complex-valued  ’’wavelet  coefficients,”  very  much 
like  the  complex  values  of  a  standard  discrete  Fourier  transform,  or  DFT.  There  are  two 
aspects  of  this  that  differ  from  standard  array- stacking  techniques,  and  which  justify  the 
extra  processing  of  the  SVD.  First,  the  SVD  looks  for  a  signal  expressed  as  a  linear  combi¬ 
nation  of  wavelet  coefficients.  Since  these  are  complex  valued,  relative  time  shifts  between 
stations  can  be  accommodated  by  the  algorithm  as  phase  shifts  in  the  wavelet  coefficients.  In 
principle,  the  SVD-decomposition  could  solve  for  the  proper  stacking  velocity  of  a  pulse,  and 
could  capture  the  non-plane-wave  nature  of  a  wave  scattered  within  or  close  to  the  sensors. 
Similarly,  the  distortion  of  the  wavefront  across  the  seismic  sensors  can  be  captured  by  the 
SVD  algorithm,  which  does  not  prescribe  a  uniform  polarization  at  all  stations.  For  Terras¬ 
cope  data,  application  of  the  SVD  algorithm  reveals  that  strong  refraction  of  surface  waves, 
usually  attributed  to  the  continent-ocean  boundary,  appears  to  occur  most  strongly  inland 
of  the  Los  Angeles  Basin  along  the  San  Andreas  Fault  (Figure  4).  Similarly,  we  observe  that 
Love-to-Rayleigh  scattered  waves  tend  to  be  stronger  for  inland,  relative  to  coastal,  Terras¬ 
cope  stations.  This  suggest  that,  from  the  point  of  view  of  forward  scattering,  the  strongest 
lateral  gradients  of  material  properties  in  the  upper  mantle  beneath  Southern  California  are 
associated  with  the  plate  boundary,  not  the  coastline. 
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Propagation  Direction  of  Rayleigh  Wave 


Propagation  Direction  of  Love  Wave 


Figure  4.  Wavelet  polarization  for  50  second  surface  waves  across  the  Terrascope  net¬ 
work  in  Southern  California.  The  large  arrow  indicates  the  great-circle  azimuth  from  the 
earthquake  source,  off  the  north  island  of  New  Zealand.  Stations  near  the  coast,  including 
SNCC,  RPV,  CALB  and  USC,  show  little  deflection  of  particle  motion,  but  large  deflections 
are  apparent  at  inland  stations  ISA,  GSC,  NEF,  VTV  and  PFO,  all  within  or  beyond  the 
San  Andreas  deformation  zone. 
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Geophys.  J.  Int.,  in  press. 

Su,  L.,  and  J.  Park,  Asymptotic  path  integral  synthetics  using  the  strong  Born  approxima¬ 
tion,  Geophys.  J.  Int.,  in  review. 

RECOMMENDATIONS  AND  FUTURE  PLANS 

By  the  end  of  the  year,  we  anticipate  submitting  manuscripts  for  publication  in  three 
areas.  In  particular,  we  plan  to  submit  manuscripts  1)  on  the  array- version  of  the  multiple- 
wavelet  transform,  comparing  crustal  surface  wave  data  from  the  Terrascope  and  Kyrghyzs- 
tan  regional  networks;  2)  on  the  hybridization  of  Love  and  Rayleigh  waves  in  1-D  anisotropic 
structures;  3)  on  the  effect  of  receiver-side  1-D  anisotropy,  with  an  arbitrary  symmetry  axis, 
on  upgoing  body  waves  e.g.  split  shear  waves  and  the  shear-coupled  PL  wave. 

Additional  areas  to  explore: 

1)  We  would  like  to  apply  the  multiple- wavelet  processing  to  an  array  with  smaller  aper¬ 
ture,  like  the  broadband  Pinon  Flat  experiment  or  the  Geyocha  array  in  Turkmenistan,  and 
to  go  beyond  a  general  study  of  wavefield  properties  to  test  the  detection  enhancement  of 
wavelet-based  techniques  on  surface  waves  at  low  signal-to-noise  ratios. 

2)  We  would  like  to  examine  the  coupling  of  anisotropic  crustal  waveguide  modes  more 
fully,  and  compare  theory  with  data  in  a  more  systematic  manner.  We  anticipate  using 
data  from  convergent-tectonic  regions,  like  the  broadband  PASSCAL  deployment  in  Tibet 
(we  are  familiar  with  this  dataset,  and  have  noted  how  Love-to- Rayleigh  scattering  becomes 
complex  as  periods  decline  from  100  to  50  seconds,  as  surface  waves  sample  the  mantle 
less  and  sample  the  crust  more.  Data  from  the  broadband  Kyrghyzstan  network  is  also 
applicable  in  this  context.  We  plan  to  test  the  behavior  at  some  US  stations  where  active- 
source  seismic  profiling  has  revealed  a  highly  reflective  lower  crust  e.g.,  within  the  Basin 
and  Range  Province.  Some  crustal  seismologists  have  attributed  this  reflectivity  to  thin 
anisotropic  layers,  which  would  develop  as  a  result  of  extensional  thinning  of  the  crust.  It 
would  be  interesting  to  see  if  coupled-mode  crustal  surface  waves  can  be  used  to  confirm  or 
deny  this  model. 
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